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ABSTRACT 

We study the buildup of magnetic fields during the formation of Population III star-forming regions, 
by conducting cosmological simulations from realistic initial conditions and varying the Jeans resolu- 
tion. To investigate this in detail, we start simulations from identical initial conditions, mandating 
16, 32 and 64 zones per Jeans length, and studied the variation in their magnetic field amplification. 
We find that, while compression results in some amplification, turbulent velocity fiuctuations driven 
by the collapse can further amplify an initially weak seed field via dynamo action, provided there is 
sufficient numerical resolution to capture vortical motions (we find this requirement to be 64 zones 
per Jeans length, slightly larger than, but consistent with previous work run with more idealized col- 
lapse scenarios). We explore saturation of amplification of the magnetic field, which could potentially 
become dynamically important in subsequent, fully-resolved calculations. We have also identified a 
relatively surprising phenomena that is purely hydrodynamic: the higher-resolved simulations possess 
substantially different characteristics, including higher infall- velocity, increased temperatures inside 
1000 AU, and decreased molecular hydrogen content in the innermost region. Furthermore, we find 
that disk formation is suppressed in higher-resolution calculations, at least at the times that we can 
follow the calculation. We discuss the effect this may have on the buildup of disks over the accretion 
history of the first clump to form as well as the potential for gravitational instabilities to develop and 
induce fragmentation. 

Subject headings: cosmology: theory — galaxies: formation — galaxies: HII regions — stars: formation 



1. INTRODUCTION 

The formation of the first stars in the Universe is a pro- 
cess well-suited to numerical computation. While direct 
observations of these stars are still an optimistic, yet- 
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unrealized prospect ([Wise fc Abel 
Catfc 



2005 



Frost et al.|2QQ9 



simulations are 

able to begin with well-posed initial conditions and di- 
rectly simulate their formation. These calculations, typi- 
cally spanning ma ny orders of magn itude in both spatial 
and density scales ( jTurk et al.|2008[|Yoshida et al.|20 08), 
include the effects of dark matter, hydrodynamics, radia- 
tive cooling, chemical heating and cooling, and the non- 
equilibrium, self-consistent evolutio n of multip le ioniza- 
tion and molecular states of the g as ([Abel et al.|199 7; AnJ 
ninoTet a l. 1997; Rip amonti fc Abeip004| IGloverpOO^ 
Glover & Abel. 2008). 

the consensus viewpoint had 
formed in isolation, between 
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l^br nearly a decade, 
been that the first stars 

30 — 300 Mq in mass, and only one per ~ 10^ M© 
halo. This had been confirmed by divergent methods: 
both Smoothed Particle Hydrodynamics (SPH) simula- 
tions and Adaptive Mesh R efinement (AMR) calcula- 
tions showed similar results (lAbel et al. 112002 lYoshida 
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aspects of this picture are being challenged. Advance- 
ments in the current generation of simulations have taken 
place in three primary areas. The first is that the effect 
of heating from the for mation of molecula r hydrogen via 
three-body reactions (iTurk et al. 2011a) is now being 
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taken into account; this will heat the gas during later 
stages of collapse and through thermal regulation of the 
accretion rate produce higher infall velocities onto the 
central core. The second aspect is related to the so- 
called "Courant myopia" of calculations proceeding to 
high densities. Because the Courant timescale becomes 
extremely short at high densities, out of necessity previ- 
ous simulations stopped after the formation of the first 
high-density core. Recent simulations have followed ac- 
cretion for up to several thousand years, using the numer- 
ical technique of sink particles, accreting Lagrangian par- 
ticles representing protostars below the resolution limit 



of the simulation (Greif et al. 2011a Clark et al. 



Stacy et al. 20101). While this approach allows to follow 



2011 



the calculations further in time it does lack the mathe- 
matical rigor that may allow one to prove the result to be 
correct. The final effect now being added is simply that 
of sampling: in the published literature, very few calcu- 
lations had been performed; compared to the wealth of 
calculations of galaxy mergers, cluster formation, and so 
on, the sampling of Population III star formation was 
dramatically underserved, with only a handful of papers 



et al.||2003||O'Shea & Norman||!^007b. However, Ti^S^f |NQrman||2007|. 



discussing even multiple simulations (see, e.g., O'Shea 



rhese recent shifts have suggested that these halos 
have fragment ed into either a small-number mu ltiplicity 



of protostars ( jTurk et al.|2009 Stacy et a l.'2010') or many 
small-mass pre-stellar objects ([CI 
|et al.||2 011a). However, " 
include! 



lark e t al. 20lf| [Greif 

none of these simulations have 
Ted the effects of magnetic fields in their calcula- 
tions, nor have they self-consistently followed the growth 
of seed magnetic fields over the collapse and virialization 
of these first minihalos. Similarly, the fragmentation has 
been only studied for times of less than 1% of the accre- 
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tion time scale of the massive stars being studied. It has 
been difficult to show whether possible early fragments 
would survive and not grow substantially in mass and 
simply merge with the central proto-star. 

The influence of magnetic fields on the collapse of the 
first stars has received a renewed interest in recent y ears 
Early numerical simulations (Xu 



jMachida et al. 



2008) 



J _ 

[et al.n200^ including the Biermann battery effect sug- 

gestea that the dynamical effect of magnetic fields on 
primordial gas is very small. Though their resolution 
was rather low, these simulations found that fields were 
primarily built up by collapse, leading to 5 ex p^l^ . The 
thermodynamics of primordial gas are quite sensitive to 
the B — p relationship v ia heating from ambipolar dif- 
fusion CSethi et al. [2010 ). Although magnetic fields are 
believed to be unlikely to affect the characteristic frag- 
mentation scale for Population III stars, they may in- 
crease the temperature in very high density g as, leading 
to higher ac cretion rates onto the protostars (Schleicher 
et al. 2009). The ability of turbulent fluid motions to 
drive dynamo action during primordial protostellar col- 
lapse has been predicted analytically b y [Schleicher et aL 



( [2010) and confirmed nuni erically by |Sur et al. ( |2Q10p 
and Federrath et al. (2011) for idealized collapse calcu- 
lations. 'I'he latter work also suggested that a minimum 
resolution of 32 elements per Jeans length is required 
to capture dynamo action, giving a very strong moti- 
vation for higher resolution cosmological MHD calcula- 
tions. Properly resolved, dynamo amplification leads to 
exponential growth of the field over a turbulent eddy 
turnover time, and can radically change the strength of 
the magnetic field at any gi ven density during collapse. 
The Federrath et al.| (|2011|) simulations considered the 
collapse of a nearly-isothermal Bonnor-Ebert sphere with 
turbulence seeded by velocity perturbations. While these 
conditions were idealized from earlier hydrodynamic cos- 
mological computations, we seek a more complete pic- 
ture by following the full magnetohydrodynamic (MHD) 
evolution of primordial gas from cosmological initial con- 
ditions. This allows us to self-consistently study the in- 
teraction of primordial gas dynamics, turbulence, and 
dynamo action. 

In this paper, we present the first highly-resolved cal- 
culations of the formation of the first stars in the Uni- 
verse from cosmological initial conditions, taking into 
account a full suite of chemical reaction rates, chemi- 
cal heating due to three-body reactions, magnetic fields 
from primordial seed fields, and a resolution of 64 zones 
per Jeans length. Furthermore, it has been conducted 
with an open source, community-built simulation code 
available for inspection and contribution. In this paper, 
we examine both the large-scale and small-scale ampli- 
fication of the magnetic field, as well as the manner in 
which resolution affects the chemo-kinetic state of the 
inner molecular cloud. A forthcoming paper will study 
the growth of turbulence in more detail. 

2. METHODS AND SIMULATIONS 

The simulations described in this work were conducted 
with the Enzo simulation cod^ Enzo is an adaptive 



^ [h ttp: //enzo-project . orgTl d escribed in [Bryan &: Nor^ 
Iman' ("lyyTfT ^U'lSJiea et ai.| ( |2UU4| ), conducted with ciiangeset 
T3cf4f 13el95 



mesh refinement (AMR) simulation code, wherein the 
baryonic fluid is discretized onto an Eulerian grid. In 
cells where certain criteria are violated (such as overden- 
sity, Jeans parameter, slopes of fluid quantities) higher 
resolution regions are inserted. This process is applied 
dynamically throughout the course of the simulation, en- 
suring adequate resolution of physical processes on all 
length scales - while the innermost regions of the calcu- 
lation may evolve more quickly dynamically, the entire 
simulation is coupled and large-scale torques, inflow, and 
halo properties are self-consistently included. 

They were initialized with cosmological initial condi- 
tions generated from Gaussian random nois e, in a man - 
ner identical to tha t conducted in Abel e ^al.| ( 2002); 
"O'Shea et ahj ( [20041 ); [Turk et al.| ( [2Q09| |20l0l ). The pa- 
rameters for generating the initial perturbations, as well 
as the partitioning of the matter budget into dark and 
baryo nic components, hav e been taken from WMAP7 re- 
sults ( Jarosik et al.|[2QTT ). The simulations were initial- 
ized from a single random seed at a redshift of 99, with 
their simulation volume centered on the first massive halo 
to form in the simulation box, 8.5 x 10^ at z = 19.72. 
For our initial conditions, we create a set of nested grids 
initialized from the primordial power spectrum, which 
are then allowed to dynamically refine as described be- 
low. The initial structure of our simulation was an out- 
ermost 128^, 300 kpc h~^ (comoving) box, within which 
three subsequent, nested levels of refinement were initial- 
ized, for a final dark matter mass resolution of 2.3M0 
and a baryonic cell width of 290 pc h~^ (comoving). 
The highest-resolution region at time of initialization is 
a cube of side length of 75 kpc h~^ (comoving) and we 
allow further refinement during the course of the simula- 
tion within a cube of side length 12 kpc h~^ (comoving) 
centered on the most massive halo in the simulation at a 
redshift of 20. 

To solve the cosmological magnetohydrodynamic equa- 
tions, we use the HLL Riemann solver with piecewise lin- 
ear reconstruction. The V • B = constraint is enforced 
using the Dedner scheme (see |Wang et al.||2008l | 2010[ 
for a full description). We have utilized a lower-oTHer 
chemical sol ver than previous studies; in both Turk et al.| 



(2010 



^ 20091 
was use^" 
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a solver based on the work of Verwer 



(1994) 



ere we have used a modified version of the 



solver from|Anninos et al. (1997|), in the interest of quan- 
tifying the influence of the magnetic field on the overall 
collapse, as well as the growth of the magnetic field. This 
choice of chemistry solver results in less finely-balanced 
equilibrium between molecular hydrogen formation and 
dissociation. Despite that, as discussed below, there are 
substantial differences in the molecular hydrogen frac- 
tion in these calculations, as well as the temperature, 
as a function of resolution. We have included the for- 
mation heating from molecular hydrogen (4.48 eV), and 
have used the Glover (2008 | three-body formation rate. 
Recent works (Naoz et aljj^lT 



Greif et al.||2 011bl) have 



suggested that relative velocities between the baryonic 
matter and dark matter may affect the details of early 
structure formation; we do not include those effects here, 
although this issue remains an important one to address 
in future simulations. 

We have conducted four simulations, each of which 
vary only in the parameters governing the AMR algo- 
rithm. Typically in AMR calculations of the first stars 
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Fig. 1. — A plot of the number of cells as a function of mass 
enclosed in those cells, for each of the simulations at the time 
at which the simulations were terminated, at a peak density of 
10~io g cm~^, the time at which all analysis in this paper was 
conducted. J 16 is represented in dotted blue, J32 in dashed green, 
J64 in solid red. This cell count is cumulative; on the left we show 
the total number of cells with mass less than the x-axis value, and 
on the right we show all cells with mass greater than the x-axis 
value. 

in the universe, Jeans resolution is the dominant refine- 
ment criteria. The four simulations presented in this 
work vary in the number of cells mandated to resolve the 
Jeans Length in each direction; we have conducted sim- 
ulations requiring 16, 32, 64 and 128 cells in each direc- 
tion, named J16, J32, J64 and J128. In all simulations, 
rather than calculate the Jeans length from the local tem- 
perature, we calculated it assuming that the ga s would 
be allowed to cool to 200 K, as in [Turk et al | ([2010j), 
which provides substantially better resolution or higher- 
temperature gas. We additionally require refinement on 
relative over densities; at the coarsest level, we fiag zones 
that are four times the mean density, and with each ad- 
ditional level this decreases by a factor of 2~^-^. This 
provides the refinement at scales larger than the virial 
radius, inside which refinement based on Jeans criteria 
dominates. These simulations began from the same ini- 
tial conditions and were seeded with an identical uniform 
magnetic field of 10"^"^ G (proper), corresponding to a 
comoving magnetic field of 10~^^ G at the start of our 
calculation. With the exception of J128, these simula- 
tions were terminated at a peak density of 10~^^ g cm~^; 
in the case of J 128, we terminated the simulation at a 
peak density of 10~^^ g cm~^. We chose this density 
as it was the point at which the simulation ceased to 
be feasible on available computational resources, and at 
which we were able to demonstrate the variation in mag- 
netic field amplification discussed below; future studies 
will utilize this level of resolution, and potentially higher, 
for calculations of Population III star formation. In Fig- 
ure [l] we plot the number of cells in each simulation as 
a function of the mass enclosed in that cell; for J64, we 
note that we have 1.4 x 10^ cells containing mass less 
than 10 times the mass of the Earth. 

3. RESOLUTION EFFECTS IN GAS DYNAMICS 

The gas dynamics of these collapsing clouds is sub- 
stantially modified by both increased resolution and the 
indirect consequences of those resolution-dependent ef- 
fects. By directly comparing the three mature calcu- 
lations (J16, J32, J64) at identical peak-density times, 
we quantify these differences. In particular, the differ- 



ences in the chemical, kinetic, and thermal states of the 
gas are the most pressing, as they will directly relate 
to the potential for fragmentation at later times. Be- 
cause these simulations are stopped before the formation 
of the first hydrostatic core, and more importantly be- 
fore that hydrostatic core has undergone several mass 
doublings, we do not use these data sets to estimate the 
initial mass function of the first stars, and instead con- 
fine our comments to measurable environmental effects 
and their potential consequences for magnetic field am- 
plification, gravitational instability, and the subsequent 
accretion onto the central molecular cloud. 

3.1. Resolution dependence of T and H2 and Vr 

In the collapse of Population III stars, the chemical 
and thermal states of the gas have been shown to be im- 
portant when considering the accretion rates and po ten- 
tial fragmentation ( Abel et al. 2002^ Q^Shea & Norman 
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et al. 



Clark et al 
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In particular, recent findings from Turk 



(2011a) suggest that the character and quantity of 



fragmentation in metal-free collapsing halos may depend 
strongly on the behavior of molecular hydrogen at high 
densities. 

As molecular hydrogen forms via three-body reactions 
(where typically the third body is either another atom of 
hydrogen or a molecule of hydrogen) the 4.48 eV binding 
energy is released into the surrounding gas in thermal en- 
ergy, via collisional deexcitation. During the collisional 
dissociation of molecular hydrogen through the inverse 
reactions to the three-body rates, the bindi ng energy is 
remov ed from the thermal energy of the gas. | Clark et al. 
( 2011 ) proposes this as a mechanism for retaining isother- 
mality in collapsing halos; as primordial gas reaches den- 
sities above 10~^ g cm~^, it becomes optically thick to 
all efficient radiative cooling until the gas reaches at least 
6000 K. This is proposed as a mechanism for allowing 
the gas to fragment under gravitational instability. This 
method of fragmentation relies heavily on the chemical 
state of the molecular cloud at high densities, and would 
be sensitive to a lower molecular hydrogen fraction, as 
this would impede gravitational instabilities from retain- 
ing isothermality and thus serving as a channel for spon- 
taneous fragmentation. 

In Figure |2] we show spherically- averaged radial pro- 
files, computed at the final output and centered at the 
densest zone in the calculation, of density, temperature, 
radial velocity, molecular hydrogen mass fraction, vor- 
ticity and enclosed mass as a function of radius. The 
central velocity was subtracted before calculation of the 
radial velocity. As a function of radius, we see good 
agreement between J16, J32 and J64 at radii greater 
than 10^^ cm in the density profile. However, at ap- 
proximately 3 X 10^^ cm (~ 1 pc) we see a small bump 
in the J 16 and J32 runs. This is not present in the J64 
runs. Typically bumps in the density plot indicate the 
formation of a second clump or other moderate level of 
fragmentation; however, based on the length and den- 
sity scales here, we attribute this to poor resolution in 
the lower-resolution runs contributing to a spike of unre- 
solved clumping just inside the virial scale; we also note 
that it correlates directly with the lowest infall velocity 
in all three simulations, so it may also simply be material 
that is being processed into the inner regions of the halo. 
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Fig. 2. — Radially-binned, spherically-averaged radial profiles of the four simulations taken at their final output. These profiles were 
centered at the most dense point; J 16 is dotted blue, J32 is dashed green, J64 is solid red, and J 128 is dot dashed cyan. Note that we 
terminated J128 before it reached the same peak density as the other simulations, although J16, J32 and J64 were terminated at roughly 
identical peak densities. Upper left is the total mass enclosed (in solar masses) at a given radius, upper right is the radial velocity in 
km s"-*^, middle left is the density inside a given radial bin in g cm""^, middle right is the temperature in K, bottom left is molecular 
hydrogen mass fraction, and bottom right is the magnitude of the local vorticity. 



The changes in the chemical and kinetic states of the 
ffas as a function of resolution are clearly correlated; with 
increased resolution, the infall velocity increases, as does 
the corresponding temperature, causing a direct decrease 
in the overall molecular hydrogen fraction. While we still 
see a peak in the molecular hydrogen, corresponding to 
a molecular hydrogen cloud of about 1 — 5 M©, the J64 
run has substantially lower molecular hydrogen fraction 
of ^ 0.2 — 0.3 in the innermost core, versus nearly twice 
that for J32 and approximately unity in the J16 case. 
As is to be expected, this correlates strongly with the 
thermal state of the gas; while J16 has a peak temper- 
ature of approximately 1400K, the J32 and J64 cases 
both exceed 2000K (with J64 reaching 2400K.) The rate 
at which molecular hydrogen is collisionally dissociated 
steeply depends on temperature, and the relationship be- 
tween^ these t wo qu antities has been explored previously 
Turk et al.| ([2011a). With increased resolution the ho- 



mogeneity of the chemical and density structure of the 
cloud cloud decreases, but we see more morphological 
homogeneity, as discussed below. 

The radial velocity structure remains roughly the same 
between all three simulations, although the magnitude of 
the infall velocity increases with increased resolution. In 
the J16 case, we see an infall velocity at the radius of the 
edge of the molecular cloud of approximately 2 km/s, 
although this increases to 6 km/s in the J64 case. The 
J32 case is in between these two. The J16 case shows 
slight peaking, indicating an under-resolved shock at the 
exterior of the molecular region, but both J32 and J64 
show very strong shocking structures. When compared 
to the enclosed mass, we see that J64 encloses a factor of 
~ 1.5 as much mass at a given radius inside the virial ra- 
dius. This likely contributes to this much greater buildup 
in shock. Furthermore, the increased resolution will al- 
low much finer shocking surfaces to develop and to be 
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Fig. 3. — The differential free-fall time for all four simulations, 
J16 (dotted blue), J32 (dashed green), J64 (solid red) and J128 
(dashed-dotted cyan, which has not reached the same peak density) 
The X-axis corresponds to a given density; the values along the y- 
axis have been calculated by subtracting the time at which the 
simulation was at a given peak density from the time at which the 
simulation had a peak density of of 10~^^ g cm~^, then dividing by 
the differential free-fall time scales between the two times. Higher- 
values therefore correspond to a longer collapse time, and lower- 
values correspond to a shorter collapse time with respect to the 
free-fall time. 

resolved (a noted problem with under-resolved hydrody- 
namics.) 

3.2. Timescales of Collapse 

Ideally, in the course of a resolution study the basis 
for comparison between two otherwise identical calcula- 
tions would be comparison at a fixed time. However, this 
method is not possible in full cosmological simulations 
because the hydrodynamics of the calculation are greatly 
influenced by resolution effects both directly, through 
turbulent stirring of the gas, viscosity, and numerical 
diffusion effects, and indirectly through variations in the 
chemical state of the e; as as a result of temperature and 
turbulence differences. 

In lieu of direct comparisons against identical times 
since the Big Bang, we have conducted our analysis 
against identical peak density times. While this neglects 
time dependencies such as the growth of mass at the 
innermost shocking region, the velocity structure varia- 
tions with time in the molecular cloud and so on, it does 
provide a baseline for comparing the chemical and kinetic 
state within the free-fallin g core. In pa r ticular , this is the 
same technique utilized in Turk et al.| (|2011a|) . 

To account for variations, we present in Figure |3] a visu- 
alization of the variations in the timescale of collapse and 
the freefall timescale for each calculation. Higher values 
correspond to a differential collapse time compared to 
expected freefall timescales; one can clearly see that sim- 
ulations J 16 and J32 are much slower to collapse than 
J64 at all densities less than 10~^^ g cm~^, at some den- 
sities even reaching ten times the freefall timescale. In 
particular, the peak delay for simulations J16 and J32 
occurs at a density of roughly 1 Q~^^ g cm~^, the density 



at which (as discussed above in § 3.1 ) molecular hydrogen 



begins forming via the three-body reaction, and therefore 



marking the edge of the molecular cloud. Furthermore, 
as noted below, the morphological differences between 
the three simulations indicate that this is correlated with 
variations in morphology. 

3.3. Morphological Differences 

In Figure [4] we show density- weighted density projec- 
tions of the final output from the calculations. From 
left, these show the J16, J32 and J64 calculations, and 
from the top they have field of view of 300 pc, 1 pc, 
and 1000 AU. The most striking difference between the 
simulations is found in the morphology. Only one (J 16) 
of the three simulations that reached 10~^ g cm~^ (J16, 
J32, J64) has formed a clear disk, with two spiral arms. 
Both J32 and J64 have not yet formed a disk, and both 
show obvious cloud-like morphology, with J64 exhibiting 
somewhat ellipsoidal morphology. While J64 will form 
an accretion disk at a later time, the variation in the set- 
tling times is striking, particularly as the time at which 
an initial protostar ignites and irradiates its surroundings 
will alter the environment for the potential formation of 
subsequent protostars. Furthermore, the dramatic differ- 
ences in the morphology suggest that resolution directly 
impacts the formation environment, and may influence 
the gravitational stability of a collapsing cloud. 

In Figure [5] we plot density- weighted projections of the 
molecular hydrogen mass fraction from the same out- 
puts and with the same field of views as in Figure [4j 
We note in particular that the morphology of the molec- 
ular hydrogen cloud roughly tracks the overall density 
structure; however, the molecular cloud in the J16 sim- 
ulation is substantially larger than those in the J32 and 
J64 calculations. This can be seen in Figure |2j wherein 
the J16 calculation has a slightly larger molecular cloud 
when averaged radially. Furthermore, while not obvi- 
ous here owing to the density- weighting, the innermost 
regions of J32 and J64 show molecular hydrogen frac- 
tions of < 0.4, and even as low as 0.2 in J64. In J16, 
the molecular hydrogen cloud tracks almost identically 
the spiral arms in density; at the leading and trailing 
ends the fraction drops. This occurs because, as the 
three-body reaction for forming molecular hydrogen is 
strongly density-dependent. The increased cooling prop- 
erties with higher molecular hydrogen fraction will only 
exacerbate this distinction between the spiral arms and 
the medium of the disk, spurring further collapse. In 
the J32 and J64 runs, we see not only highly-irregular, 
spheroidal structure, but we see little to no evidence for 
runaway gravitational instability at this time in the col- 
lapse; in fact, the roughly spheroidal structure of the 
cloud suggests that disk fragmentation at the '^100 AU 
scale may be disfavored, certainly until a later time. 

4. MAGNETIC FIELD AMPLIFICATION 

We infer magnetic field amplification above that aris- 
ing simply from the spherical compression of frozen-in 
field lines during collapse by plotting the magnetic en- 
ergy, Eb = B'^/Stt, as a function of density for each of 
our four simulations (figure [6|. Assuming a power law 
relation between magnetic energy and density, Eb oc p^, 
a spherical collapse would result in b = 4/3. Any steeper 
values of b must result in additional amplification in the 
form of a dynamo. We expect a small-scale turbulent 
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Fig. 4. — Density-weighted projections (through the entire simulation domain) of the average density field in simulations J16 (left 
column), J32 (second column) and J64 (third column) at fields of view of 300 pc (top row), 1 pc (middle row), and 1000 AU (bottom row). 



dynamo if there is significant kinetic energy in turbulent 
motions. Figure |6] shows that for the best resolved sim- 
ulations, J64 and J128, the magnetic field scales with a 
much higher power, 6 1.78 at all densiti es p > 1Q~^^, 
roughly consistent with the findings of |Federrath et al.| 

1 2Q11| ). By contrast, the lower resolution runs (J 16 and 
32) begin to show 6 > 4/3 behavior at the same den- 
sity as the higher resolution runs, but never reach the 
same b as the higher runs and eventually begin to show 
a shallower slope of Eb with p at the highest densities 
{p > 10~^^). In the outer part of the collapse, near the 
virial radius, all runs show b 4/3, suggesting that at 
low densities, the growth of magnetic energy is primarily 
due to the roughly spherical collapse, though the data is 
somewhat scant. 

4.1. Velocity Structure In the Collapse Region 



The dynamo activity in the previous section is reso- 
lution dependent, but robust above roughly 64 cells per 
Jeans length (J64); that is, we continue to see more field 
amplification, consistent with a minimum resolution re- 
quirement of between 32 and 64 zones per Jeans length. 
This is broadly consonant with the results of Federrath 
et al, who find a similar resolution cutoff in their simu- 
lations of magnetic field growth in the nearly-isothermal 
collapse of a Bonnor-Ebert sphere. In their picture, the 
reason for this cutoff is the lack of sufficient power in ro- 
tational motions when the Jeans length is resolved with 
fewer than ~ 30 cells. Our resolution requirement for 
strong and sustained dynamo action appears to be about 
a factor of two larger, which could be due to our use of 
a considerably more diffusive MHD solver (we use the 
HLL solver; they use HLL3R). In order to solidify the 
connection between dynamo action and the presence or 
absence of turbulent fluctuations, we flrst consider pro- 
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Fig. 5. — Density- weighted projections (through the entire simulation domain) of the average molecular hydrogen mass fraction field in 
simulations J 16 (left column), J32 (second column) and J64 (third column) column) at fields of view of 300 pc (top row), 1 pc (middle 
row), and 1000 AU (bottom row). 

magnetic energy growth with density begin to fah off in 
low resolution simulations. This correlation between vor- 
ticity and magnetic energy is consistent with small-scale 
dynamo action generated by incoherent velocity fields. 



jections of uj'^ where = V x v is the fiuid vorticity cen- 
tered at the densest point in the computation (figure [s]). 
At the 300pc scale (figure [sj top row), the vorticity is 
identical among the four resolutions. However, zoom- 
ing in to ~ Ipc, there are already significant differences 
between all four runs. For the other three, a trend is 
clear: increasing resolution creates much larger regions 
of high vorticity that is indicative both of an overall in- 
crease in the turbulent energy but also of a d ecreasing 
coherence of the collapse region (see section 3.3). Finally, 
at lOOOAU, the vorticity structure is completely different 
between the ordered disk-like structure in J16 and J32 
and the amorphous turbulent core of J64. It is thus clear 
that as a function of resolution, we see an increase in 
vorticity production and a decrease in the characteristic 
scale of that vorticity at the smallest scales within the 
collapsing core. These high density regions are just where 



4.2. Magnetic Field Saturation 

The kinematic, small-scale dynamo acts by twisting 
magnetic field in a turbulent fiow field. The random- walk 
character of the fiow will lead to a continuous stretching 
and thus strengthening of the magnetic field. In a typical 
numerical model of the small-scale dynamo, turbulence 
is driven in a fiuid with a tiny seed field; the field grows 
exponentially until it nears equipartition with the turbu- 
lent velocity field, at which point the Lorentz force reacts 
back on the fiuid, leading to a strongly nonlinear coupling 
between v and B which ultimately saturates the growth 
of the magnetic field. 

Here, the collapse timescale is decreasing with time. 
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p (g cm ) 

Fig. 6. — Lines show volume- weighted magnetic energy scaled by 
p^l^ as a function of density for a 1 — kpc sphere centered on the 
densest point in the final output for each of four resolutions, J16, 
J32, J64, and J128. The colored pixels give the mass in Eb — p bins 
for the J64 run. We have overlaid two power laws, giving B oc p^/^ 
(lower dashed line), the expected scaling for spherical collapse of 
frozen-in magnetic fields, and B oc p-*^-^^ (upper dashed line), the 
best-fit power law to the J64 data. 

which means that the computation reaches its end long 
before the magnetic field reaches saturation. We stop 
each of our runs when their peak density is pmax — 2.2 x 
10-^°, except for the J128 run, which was terminated 
at pmax — 2.0 X 10~^^ due to constraints on computing 
time. Nevertheless, even in the J64 run, where dynamo 
action is most vigorous, the kinetic energy exceeds the 
magnetic energy in mass weighted averages by a factor 
of - 3600. 

The growth rate a of the small-scale dynamo is a func- 
tion of the Reynolds number Re, the ratio of the advec- 
tive timescale to the viscous timescale, which we do not 
control in these runs (dissip ation is solely a nume rical 
artifact of finite resolution). ( [Federrath et al.||2011 see, 
e.g.,) found the growth rate to be a cx Re " . Even if 
we were to solve the Navier- Stokes equation including 
viscosity and the non-ideal induction equation includ- 
ing Ohmic resistivity, we would be nowhere near the ac- 
tual Reynolds numbers of the fiows in primordial gas. A 
rough estimate of Re using the Spitzer viscosity, charac- 
teristic velocity and length scales taken from Figure [2] is 
about Re ~ 10^^ at a radius of ~ 1000 AU . This means 
that our simulations cannot make a definitive statement 
regarding the saturation level of the collapse dynamo, 
though our results can be taken as a lower limit to the 
field strength during collapse. We note that this indi- 
cates that magnetic fields may be dynamically important 
in the formation of the first stars, and that their amplifi- 
cation requires much higher resolution than is currently 
being utilized. 

5. DISCUSSION 

It is clear from our simulations that gravitational col- 
lapse from cosmological initial conditions does indeed 
drive a small-scale dynamo if the fiow is resolved well 
enough to capture the turbulent fiuct nations. This tur- 
bulence is produced self-consistently by the collapse dy- 
namics, and is not seeded by any driving or initial per- 



turbations. The change in slope between J64 and J128 
indicates that while we have demonstrated the minimum 
resolution for dynamo action to occur, further amplifi- 
cation of the magnetic field will result from increased 
resolution, and thus our results should not be considered 
converged. We will explore the details of the generation 
the turbulence and the details of the dynamo generated 
magnetic field in a forthcoming paper. 

The most striking difference between the J 16, J32 and 
J64 runs is the resolution-dependence in the chemical and 
kinetic states of the gas. While we initially expected that 
the simulations would produce largely the same results 
in the structure of the collapsing clouds, the three sim- 
ulations produced substantially different disk structure. 
Because we examined these three simulations at iden- 
tical peak density values, the time remaining until the 
formation of the first hydrostatic core should be roughly 
identical in each run; as such, the initial structure of 
th e cloud will like ly not change in that intervening time. 
, Clark et al. (2011 ) presented simple comparisons between 
the accretion rate onto the central accretion disk and the 
accretion rate onto a central protostar; in those calcula- 
tions these comparisons indicate that the central core will 
be unable to process sufficient gas to deplete the disk, 
which they observed result in gravitational instability. 
We note that in our simulations, we see substantially re- 
duced disk-like structure, as well as substantially higher 
infall velocities at the edge of the molecular hydrogen 
cloud. The combination of these two effects results in an 
unclear change to the potential fragmentation character- 
istics of these collapsing clouds; reduced angular momen- 
tum should allow greater accretion onto a central core, 
but the increased radial velocity onto the molecular cloud 
will result in greater mass-buildup. Furthermore, we note 
that while we do not resolve the formation of a central 
core or the subsequent accretion, and thus cannot make 
statements regarding the fragmentatio n characteristics in 
our simulation, the result presented in Clark et al. (2011 ) 
was constructed from a series of fixed resolution calcula- 
tions with a much lower Jeans resolution than in our J64 
run (and indeed our J32 run); the results in this paper 
suggest that higher resolution may change the specific 
nature of fragmentation. 

The correlation between the dissociation of molecular 
hydrogen and the increased resolution was unexpected, 
but in retrospect is not entirely surprising. With in- 
creased infall velocity, increased turbulent support, and 
increased overall temperature, the rate at which molec- 
ular hydrogen collisionally dissociates increases steeply, 
while the rate at which it associates via three-body re- 
actions decreases. With increased turbulent energy, and 
a corresponding increase in the temperature, the inner- 
most region in the pre-stellar molecular cloud will natu- 
rally undergo some dissociation. However, what remains 
unclear is how this will affect the response o f the gas to 
gravit a tional instability. P rior to the work of jClark et al. 
(2011); Greif et al. (2011a ), fragmentation had been seen 
at scales of ^ 2000 AU ( [Turk et al.||2"QQ9| IStacy et al.] 
2010); at smaller scales, fragmentation had (prior to the 
works of Clark, Greif and their collaborators) been be- 
lieved to be suppressed because of the inefficient cool- 
ing of molecular hydrogen. In these more recent works, 
however, the ability of the gas to shed thermal energy 
through dissociation of molecular hydrogen, essentially 
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Fig. 7. — Density-weighted projections (through the entire simulation domain) of the average magnetic energy in simulations J16 (left 
column), J32 (second column) and J64 (third column) at fields of view of 300 pc (top row), 1 pc (middle row), and 1000 AU (bottom row). 



allowing compressional heating to be freely reversed, has 
been identified as the mechanism by which fragmenta- 
tion proceeds. At high densities (the densities at which 
sink particles have been inserted in their calculations), 
the gas becomes effectively isothermal (and later, once 
the molecular hydrogen has been depleted, adiabatic.) 
However, this hypothesis relies on a substantial reservoir 
of molecular hydrogen; while a gas cloud containing a de- 
pleted but non-zero supply of molecular hydrogen is still 
likely to undergo a quasi-isothermal phase during its col- 
lapse, the resilience of that gas and the stiffness of its 
effective equation of state are likely to be modified. Fur- 
thermore, because the fraction of molecular hydrogen, 
and thus the equation of state of potentially-fragmenting 
gas sensitively depends on the rate at which molecular 



hydrogen is fo rmed and destroyed, as discussed in [Turk 
[et al.| ( (2Qlla[ ), this provides an additional uncertainty. 
Without extremely high-resolution calculations that fol- 



low the formation of the first core, as well as its subse- 
quent accretion over several thousand years, the charac- 
ter of fragmentation in accretion disks around Population 
III pre-stellar cores is still an open question. 

6. CONCLUSIONS 

We have presented on the first fully-cosmological calcu- 
lations of Population III star formation that include all 
relevant chemical processe s, as well as magnetic fields. 



We see, in agreement with Federrath et al. (2011), that 



a critical resolution exists above which we are able to 
resolve small-scale dynamo action resulting in increased 
magnetic-field field growth. While we adequately resolve 
the action of this small-scale dynamo, we have not yet 
resolved the amplification caused by that dynamo, as 
demonstrated by the nascent J128 simulation. In fact, 
while these calculations have yet to show a dynamically- 
important magnetic field, estimates of the possible sat- 
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uration level of the magnetic field indicate that with in- 
creased resolution or stronger seed fields, magnetic fields 
may in fact become dynamically important. 

A side effect of conducting these resolution studies has 
been that we see substantial variation in the chemical, 
kinetic, and velocity structure of collapsing metal-free, 
star- forming clouds. In particular, the substantial dif- 
ference in the J64 run suggests that under-resolving the 
hydrodynamics can result in incorrect solutions: the frac- 
tion of hydrogen gas that is in molecules, the speed of 
sound, the infall velocity and the turbulent support all 
depend strongly on hydrodynamic resolution. 

Our calculations, while falling short of following po- 
tential fragmentation in detail, suggest that resolution 
effects may be an additional complication in determin- 
ing the initial mass function of the first stars; in fact, 
they suggest that previous calculations in the literature 



have not yet converged on a solution. Future calcula- 
tions, of much higher resolution and attaining a greater 
saturation level of magnetic fields, will help to resolve 
these outstanding questions. 
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